library(tidyverse)
library(ggpubr)

pub_sle <- read_delim('~/append-ssd/nextflowing/sle_vcf/imputed/pub_sle_imputed.vcf.gz', comment = '##')
pub_panel_sle <- read_delim('~/append-ssd/nextflowing/sle_vcf/imputed/pub_sle_paneled_imputed.vcf.gz', comment = '##')

plink_filter <- read_delim('~/append-ssd/nextflowing/sle_vcf/imputed/plink.export.vcf.gz', comment = '##')

plpub_sle |> filter(POS == 105737776) -> igg1_sle

igg1_sle <- igg1_sle |>
  select(c(POS, INFO, CTRL_10:last_col())) |>
  pivot_longer(CTRL_10:last_col(), names_to = 'sample', values_to = 'genotype')

igg1_sle <- igg1_sle |>
  mutate(genotype = str_remove(genotype, '.+:') |> as.numeric())

igg1_sle <- igg1_sle |>
  mutate(group = case_when(str_detect(sample, 'MC|FC|healthy|HC|CTRL') ~ 'ctrl', .default = 'SLE') )

igg1_sle |> ggplot(aes(group, genotype)) + geom_violin() + stat_compare_means()
igg1_sle |> summarise(mean(genotype), .by = group)

igg1_sle |>
  mutate(`#IID` = sample, SLE = ifelse(group == 'SLE', 2, 1), .keep = 'none') |>
  relocate(`#IID`) |>
  write_tsv('pub_sle.plink2_pheno')

# pca results of plink2
plink.pca <- read_delim('~/append-ssd/nextflowing/sle_vcf/imputed/plink2.eigenvec')

plink.pca |> left_join(igg1_sle, join_by(`#IID` == sample)) |>
  ggplot(aes(PC1, PC2, color = group)) + geom_point()

# glm results of plink2
plink.glm <- read_delim('~/append-ssd/nextflowing/sle_vcf/imputed/plink2.SLE.glm.logistic.hybrid') |>
  filter(ERRCODE == '.')

plink.glm |>
  filter(POS == 105737776)

plink.glm |>
  filter(P < .05)
